# Indicates which parties are still alive in the next election at t + 1. 
partydied <- function(x) {
	elections <- unique(x$year)
	n <- length(elections)
	for(i in 1:n){
		if (i != n){
			now.part <- x$party[x$year==elections[i]]
			tom.part <- x$party[x$year==elections[i+1]]
			x$die[x$year==elections[i]][!(now.part %in% tom.part)] <- 1
			}
		}
	return(x)
	}


# Calculates lost/gains in seats. 
seatsdiff <- function(x) { 
	if (nrow(x) > 1){
		x <- x[order(x$year), ]
		n <- length(x$seats)
		x$seatsdiff <- (x$seats - c(NA, x$seats[1:(n-1)]))
		}
	return(x)
	}

# Calculates which party received most gains in an election. If more than 
# 50% of the parties compete for the first time in an election, the winner 
# of the election is choosen instead of the party with the most gains. 
maxgain <- function(x) {
	if ( sum(is.na(x$seatsdiff)) < 0.5 * nrow(x) ){
		x$gain[x$die==0][which.max(x$seatsdiff[x$die==0])] <- 1
		}
	else{ x$gain[which.max(x$seats)] <- 1	}
	return(x)
	}

# Builds election bridge (incentive hypothesis).
electionbridge1 <- function(x){
	x <- x[order(x$year),]
	x$t <- 1
	if ( nrow(x) != 1) {
		for(i in 2:nrow(x)){
			x$t[i] <- x$t[i-1] + ifelse( ( x$gain[i] ==0 & x$gain[i-1] ==0) | ( x$gain[i] ==1 & x$gain[i-1] ==0), 1, 0)
			}
		}
	x$T <- max(x$t)	
	return(x)
	}

# Builds election bridge (no hypothesis).
electionbridge2 <- function(x){
	x <- x[order(x$year),]
	x$t <- 1
	if ( nrow(x) != 1) {
		for(i in 2:nrow(x)){
			x$t[i] <- x$t[i-1] + 1
			}
		}
	x$T <- max(x$t)	
	return(x)
	}

# Builds election bridge (winner hypothesis).
electionbridge3 <- function(x){
	x <- x[order(x$year),]
	x$t <- 1
	if ( nrow(x) != 1) {
		for(i in 2:nrow(x)){
			x$t[i] <- x$t[i-1] + ifelse( ( x$won[i] ==0 & x$won[i-1] ==0) | ( x$won[i] ==1 & x$won[i-1] ==0), 1, 0)
			}
		}
	x$T <- max(x$t)	
	return(x)
	}

# Builds cross-national bridge.
nationalbridge <- function(x){
	if(nrow(x) > 1){
		x[x$study==2, c("t","T", "i")] <- x[x$study==1, c("t","T", "i")]
		}
	else{x$skp <- 1}
	return(x)
	}


# Builds EP party intercept indicator. 
b0indicator <- function(x){
	if(1 %in% x$overlap){ 
		x$z[x$overlap==0] <- as.numeric(x$it[x$overlap==1]) 
		}
	else{ return(x) }
	return(x)
	}

# Makes an index variable named 'to' from the variables in 'from'.
indexmaker <- function(data, from, to){
	if( length(from) == 2 ){ data$ind <- as.factor(paste(data[, from[1] ], data[, from[2] ], sep="_")) }
	if( length(from) == 1 ){ data$ind <- as.factor(paste(data[, from ], sep="_")) }
	if( length(from) > 2 ){ stop("'from' has to be of length 2.") }
	levels(data$ind) <- seq(1, nlevels(data$ind))
	colnames(data) <- c(colnames(data)[1:(ncol(data)-1)], to)
	return(data)
	}

# Re-indexing a variable. 
reindexing <- function(data, which){
	data[,which] <- as.factor(as.numeric(data[,which]))
	levels(data[,which]) <- seq(1, nlevels(data[, which]))
	return(data)
	}
	
# Builds parameter selection datasets. 
paramselect <- function(data){
	
	phi.select <- subset(data, select=c("countryname", "entity", "year", "party" ,"it","i","t","T", "w","z"))
	phi.select <- phi.select[duplicated(phi.select[,c("it","i","t","T")])==FALSE,]
	phi.select <- phi.select[order(phi.select$it), ]
	phi.select$z[phi.select$z==0] <- phi.select$w[phi.select$z==0] + nrow(phi.select)

	theta.select <- subset(data, select=c("entity", "year", "ce"))
	theta.select <- theta.select[duplicated(theta.select)==FALSE, ]
	theta.select$start <- 0
	theta.select <- ddply(theta.select, "entity", function(x) {x$start[which.min(x$year)] <- 1; return(x)} )
	theta.select <- theta.select[order(theta.select$ce),]

	rho.select <- subset(data, select=c("entity", "c"))
	rho.select <- rho.select[duplicated(rho.select)==FALSE, ]
	rho.select <- rho.select[order(rho.select$c),]
 	
	out <- list(phi.select=phi.select, theta.select=theta.select, rho.select=rho.select)
	return(out)
	
	}

# Gives back the cross-national bridges. 
whereisnationalbridge <- function(x){ c( length(unique(x$party)), unique(x$year[x$study==1]), unique(x$year[x$study==2])) }

# Calculates the range of a variable.
rangedist <- function(x){
	min <- range(x, na.rm=TRUE)[1]
	max <- range(x, na.rm=TRUE)[2]
	dist <- max - min
	return(dist)
	}

# Indicates if x[1] and [2] overlap 0.
overlap0 <- function(x)	{
	x <- sort(x)
	is <- ( x[1] < 0 & x[2] > 0  )
	return(is)
	}

# Extracts a MCMC sample. 
getMCMCSample <- function(mcmclist,i){
  chainext <- function(x,i) return(x[,i])
  return(as.mcmc.list(lapply(mcmclist, chainext, i = i)))
  }

# Runs Gelam-Rubin on an MCMClist.  
makeGelmanDiag <- function(mcmclist){
  n <- dim(mcmclist[[1]])[2]
  gd <- matrix(NA,n,2)
  for(i in 1:n){
    mcmcsample <- getMCMCSample(mcmclist,i)
    gd[i,] <- t(as.vector(gelman.diag(mcmcsample)$psrf))
    }
  colnames(gd) <- c("Point Est.", "Upper C.I")
  gd <- data.frame(parm=colnames(mcmclist[[1]]), gd)
  return(gd)
  }

# Calculates phi's mean / BCI from a model posterior. 
doresults <- function(x, datafolder=NULL, ov=FALSE){

	if( is.null(datafolder) ) datafolder <- x

	# Load posterior. 
	load(paste("./models/",x,"/out_jags.Rdata",sep=""))

	# Remove burn-in.
	phimcmc <- mcmc.list(mcmc(out_jags$phimcmc[[1]], 1001,2000), mcmc(out_jags$phimcmc[[2]], 1001,2000))
	phimcmc <- as.matrix(phimcmc[,str_detect(colnames(as.matrix(phimcmc)), "phi"),])

	# Delete prior values.
	k <- dim(phimcmc)[2] - 8
	phimcmc <- phimcmc[,1:k]	

	# Calculate mean / quantiles. 
	phimu <- apply(as.matrix(phimcmc), 2, mean)
	phiq <- apply(as.matrix(phimcmc), 2, quantile, probs=c(0.025, 0.975))

	# Calculate overlap with prior distribution
	if ( ov ==TRUE ){

		load("chesprior.Rdata")
		sigmainv0_ <- chesprior$sigmainv0[out_jags$phi.select$w] 
		mu0_ <- chesprior$mu0[out_jags$phi.select$w]

		index <- as.matrix(seq(1,length(mu0_)))
		phisim <- t( replicate( nrow(phimcmc), apply(index, 1, function(x) rnorm(1, mu0_[x], sqrt(1/sigmainv0_[x]) ))) )

		diff <- (phimcmc - phisim)
		diffq <- apply(diff,2, quantile, probs=c(0.05, 0.95))	
		ovl <- as.numeric(apply(diffq, 2, overlap0))

		phisum <- cbind(out_jags$phi.select[,c("it")], phimu, t(phiq), ovl)
		colnames(phisum) <- c("it", paste("lr.mu.",x,sep=""), 
						  				  	 paste("lr.lo.",x,sep=""), 
											 paste("lr.up.",x,sep=""),
											 paste("lr.ov.",x,sep=""))
	} else {
		phisum <- cbind(out_jags$phi.select[,c("it")], phimu, t(phiq))
		colnames(phisum) <- c("it", paste("lr.mu.",x,sep=""), 
						  				  	 paste("lr.lo.",x,sep=""), 
											 paste("lr.up.",x,sep="") )
	}
	# Merge with original data.												
	load(paste("./models/", datafolder, "/for_jags.Rdata", sep=""))	
	cmpemp <- for_jags$cmpemp
	phisum <- merge(cmpemp[,c("countryname", "party", "year", "study", "it")], phisum, by="it")
	phisum$it <- NULL

	# Save.
	save(phisum, file=paste("./models/", x, "/meanbci_phi.Rdata", sep=""))
	}
	
	
# Runs Gelman-Rubin on the output of a model.
dodiags <- function(x){
	file <- x
	load(paste("./models",file,".Rdata",sep=""))
	phimcmc <- out_jags$phimcmc
	auxmcmc <- out_jags$auxmcmc

	gd_auxmcmc <- makeGelmanDiag(auxmcmc)
	gd_phimcmc<- makeGelmanDiag(phimcmc)

	cat("\n")
	cat(  "Gelman-Rubin Diagnostics", append=F, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))
	cat(  "========================", append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))

	cat(  "\nphi \n", append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))
	out <- capture.output(summary(gd_phimcmc))
	cat(  out, append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))

	cat(  "\nb \n", append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))
	out <- capture.output(summary(gd_auxmcmc[str_detect(gd_auxmcmc$parm, "^b"), ]))
	cat(  out, append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))

	cat(  "\nsigma\n", append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))
	out <- capture.output(summary(gd_auxmcmc[str_detect(gd_auxmcmc$parm, "^sigma"), ]))
	cat(  out, append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))

	cat(  "\nlambda", append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))
	out <- capture.output(summary(gd_auxmcmc[str_detect(gd_auxmcmc$parm, "^lambda"), ]))
	cat(  out, append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))

	cat(  "\ntheta\n", append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))
	out <- capture.output(summary(gd_auxmcmc[str_detect(gd_auxmcmc$parm, "^theta"), ]))
	cat(  out, append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))

	cat(  "\nrho\n", append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))
	out <- capture.output(summary(gd_auxmcmc[str_detect(gd_auxmcmc$parm, "^rho"), ]))
	cat(  out, append=T, sep="\n", file=paste("./models",file,"_diags.txt", sep=""))
	}
